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Abstract 

We identify ghost trajectories of symplectic-energy-momentum (SEM) 
integration and show that the ghost trajectories are not time reversible. 
We explain how SEM integration can be regularized, in a SEM pre- 
serving manner, so that it is time reversible. We describe an algorithm 
for implementing the regularized SEM integrator. Simulation results 
for the pendulum are given. Coordinate invariance of the regularized 
SEM integrator is briefly discussed. 

Key Words DTH dynamics, symplectic energy momentum integra- 
tor, discrete mechanics, discrete time Hamiltonian, discrete vari- 
ational principles, principle of least action, energy conserving 
methods, extended phase space, midpoint method. 



1 Introduction 

Is symplectic-energy-momentum (SEM) integration obstructed by singular- 
ities? In 1994, the answer appeared to be yes. SEM integration of the rota- 
tions of a pendulum did not appear possible, Shibberu |18j . In this article, 
we show that SEM integration is actually not obstructed by singularities. 
Difficult to compute, ghost trajectories are identified, but these ghost tra- 
jectories are not time reversible. We explain how SEM integration can be 
regularized, in a SEM conserving manner, so that it is time reversible. We 
begin with a brief review of SEM integration. 
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A SEM integrator is a symplectic integrator which exactly conserves 
energy and momentum. Symplecticness implies the integrator can be derived 
from a discrete variational principle |2L)| . The discrete variational principle 
of a symplectic integrator gives it coordinate invariant properties. 

SEM integration emerged from two lines of research, symplectic inte- 
gration and discrete mechanics. Efficient computation is emphasized in 
symplectic integration. Preservation of the physical laws of nature is the 
emphasis in discrete mechanics. The term "symplectic-energy-momentum 
integrator" was coined and popularized by Kane, Marsden and Ortiz |12j . 
See also Chen, Guo and Wu [3] for related work on higher-order, sympletic- 
energy integrators. Guibout and Bloch ^Tj have developed a general frame- 
work for deriving many of the published symplectic integrators, including 
SEM integrators. They also provide an interesting comparison of the various 
discrete variational principles used. 

The author's work on a discrete-time theory for Hamiltonian dynamics 
(DTH dynamics) predates the work of Kane et. al. [T^]. DTH dynamics 
is a SEM integrator. DTH dynamics originated from an effort to obtain the 
exact energy and momentum conserving properties of the discrete mechanics 
of Greenspan [S], ^U] from the variational principle used in the discrete 
mechanics of Lee [Hj, |15j . DTH dynamics was proved in 1994 (see Shibberu 
|18j . |19j ) to be symplectic and hence a SEM integrator. D'Innocenzo, Renna 
and Rotelli |Sj have done work that can also be related to SEM integration. 

In the author's work, SEM integration is accomplished by varing the 
time step of the midpoint scheme to enforce exact energy conservation at 
the midpoint of each step. Symplecticness and momentum conservation 
occur at the vertices of each step. The relationship between the time step 
and energy conservation originates from the fact that the negative of the 
energy (Hamiltonian) is the momentum corresponding to time, Shibberu 

na, Lee m. 

The requirements of symplectic-energy integration are highly restrictive 
as illustrated by Ge's Theorem jjj, j^j. An early existence and uniqueness 
result for DTH dynamics was given in Shibberu and an explanation of 
why Ge's Theorem is not violated was given in Shibberu JH]- The suffi- 
cient condition for the existence and uniqueness of DTH trajectories proved 
in Shibberu ^Zj does not cover all the points in phase space where the 
Hamiltonian function is smoothly defined. It first appeared, from simula- 
tion results in Shibberu ^E]; that DTH dynamics was obstructed by points 
where this sufficient condition does not hold. Kane et. al. [12] later observed 
similar difficulties near points they refer to as "turning points" and Chen et. 
al. |3] also mentioned the need to avoid singularities in their algorithm. In 
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this article, we explain how SEM integration can be regularized in a manner 
which preserves SEM properties and time reversibility at such points. 

The outline of this paper is as follows. In section |5J we review the foun- 
dations of DTH dynamics. We introduce a discretization of Hamiltonian 
dynamics which is equivalent to, but simpler than the discretization used 
in Shibberu An example of a ghost trajectory of the pendulum is 

also given. In section 03 we introduce the two complementary variational 
principles used to regularized DTH dynamics. Regularized DTH equations 
are derived and shown to preserve SEM properties. Coordinate invariance 
of the regularized equations is briefly discussed. In section |1J we give a 
detailed description of an algorithm for solving the regularized DTH equa- 
tions. Numerical results for the pendulum and Kepler's one-body problem 
are discussed. Finally, certain peculiarities of regularized DTH dynamics 
are described. 



2.1 Foundations of DTH Dynamics 

We begin by introducing an extended-phase space version of the principle of 
least action. Let H(t, q±, . . . , q n , pi, . . . ,p n ) be the Hamiltonian function of 
an n-degree of freedom Hamiltonian dynamical system where t is time and 
qi and pi, i = 1, ... ,n, are position and momentum coordinates. Let q = 
(qi, . . . , q n , t) T andp = (p\, . . . ,p n , p) T be extended phase space coordinates 
where p is the momentum conjugate to time. (See (TH] , [H] or JH] for a 
description of p.) Let z = (q,p) J ■ Consider the extended-phase space action 
integral 



and / is the n + 1 dimensional identity matrix. The extended-phase space 
Hamiltonian function is H(z) = p+H(t, q±, . . . , q n ,Pi, ■ ■ ■ ,Pn)- The principle 
of least (stationary) action states that the trajectory z(t) of a Hamiltonian 
dynamical system cause the action integral A(z(-)) to be stationary under 
variations which satisfy the boundary constraints q(ro) = qo, p(rf) = pf and 
the Hamiltonian constraint 7~t(z) = 0. (Additional details are given in |19j.) 

The action integral A(z(-)) is discretized in 19 a by evaluating the integral 
along piecewise-linear, continuous trajectories and then appending boundary 
terms to account for the boundary conditions. An equivalent discretization 
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Figure 1: A piecewise- linear, continuous trajectory in extended-phase space. 

with no boundary terms is described below. This discretization makes it 
possible to provide a simpler derivation of the DTH equations. 

Lemma 1 Let da k represent the boundary of the triangle a k in extended- 
phase space shown in Figure^ Then, along dak, 



Lemma n follows from Stoke's formula pQ. 

Definition 2 (One Step Action) The one-step action of a discrete-time 
Hamiltonian dynamical system is defined to be 



A discrete-time Hamiltonian (DTH) trajectory is a piecewise-linear, con- 
tinuous trajectory which satisfies the following discrete variational principle. 

Definition 3 (DTH Principle of Stationary Action) The one-step ac- 
tion A(zk, Zk+i) is stationary along a DTH trajectory for variations which 
fix qk and Pk+i and satisfy the Hamiltonian constraint ri(zk) = where 
Zk = \{zk+i + Zk) and k = 0, . . . , N - 1. 

Theorem 4 (DTH Equations) A DTH trajectory is determined by the 
following equations: 



Az(-)) = [ 1 -z(t) t Jz'{t) t dr = lAq k T A P k 



Jdcr k 



A{zk,z k+ i) = -Aq k Ap, 



> k , where z k = (q k ,Pk) 



Az k = X k JH z (z k ) 
H(z k ) = 0. 



(la) 
(lb) 
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The proof of Theorem 0] follows from the proof of Theorem |f)J The proof 
that equations (|laj l — l|lb |) is a SEM integrator (i.e. preserve symplecticness 
and conserve momentum at Zk and conserves energy at Zk ) can be found in 



2.2 An Example of SEM Integration 

The extended-phase space Hamiltonian function for a pendulum is TC(q,p, p) 



p + 2P ~ cos(q). The corresponding DTH equations are 

Aq k = X k Pk 

Ap k = -A fc sin(c/ fc ) 
Apfc = 

= Pk + TjPk ~ cos (Qk) 

Observe from equation ()2b|) that equals the time step Atk- 



(2a) 
(2b) 
(2c) 
(2d) 

(2e) 



Figure 2(a) shows two DTH trajectories projected onto the phase por- 
trait of the pendulum. Observe that the linear segments of DTH trajectories 
are tangent to their respective energy-conserving manifolds (except possibly 
where they cross the v-shaped curves). 

A sufficient condition for the existence and (local) uniqueness of so- 
lutions to equations (|2*a|) - ((2*e| ) is the condition ip(z) 7^ where tp(z) = 
( JTl z ) T 7i zz ( J7i z ) El- The function tp(z) plays a key role in the regulariza- 
tion of SEM integration. The v-shaped curves in Figure [2(a) are points in 
phase space where ip(z) =0 and were first described in |18j . 



2.3 Ghost Trajectories 



Figure 2(b) shows that the DTH equations (|2aj) -(|2e | ) are poorly conditioned 
when tp is near zero. Never the less, these equations determine trajectories 



which cross the v-shaped curves in Figure 2(a) in a SEM conserving manner. 
However, these trajectories are not time reversible as is seen in Figure|21 The 
linear segment crossing ip(z) =0 forward in time is tangent to the energy- 
conserving manifold on the left side, but the segment crossing backward 
in time is tangent on the opposite side. We will call these non-reversible 
trajectories, ghost trajectories. The trajectory shown crossing the v-shaped 



curves in Figure 2(a) has been regularized in a SEM conserving manner so 
that it is time reversible. How this is done is explained in the next section. 
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(a) Phase portrait 



(b) Condition number 



Figure 2: DTH trajectories for the nonlinear pendulum. The v-shaped 
curves correspond to points where ip = 0. 
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(a) View of H = 0, tp = 



(b) Magnified view of H = 0, 
ip = 



Figure 3: A ghost DTH trajectory crossing ip = forward in time (solid 
curve) and then time-reversed so that it crosses backward in time (dashed 
curve) . 
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3 Regularized DTH Dynamics 



3.1 Regularization 

We regularize DTH dynamics by resorting to two complementary variational 
principles. The first variational principle restricts variations in Definition 
121 by the inequality constraint ip(z k ) > ipk where ip k is a constant. The 
second variational principle restricted variations by the inequality constraint 
ip(zk) < ipk- We alternate between the two variational principles to generate 
trajectories which cross ip = in a time reversible manner. 

Definition 5 (Regularized DTH Principle) The one-step action A(z k , z k +i) 
is stationary along a DTH trajectory for variations which fix q k andp k+ i and 
satisfy the Hamiltonian constraint TL(z k ) = and the inequality constraint 
ip(zk) > ipk (or ip(z k ) < ip k ). 

Theorem 6 (Regularized DTH Equations) Assume TL z (z k ) andip z (z k ) 
are linearly independent for k = 0, . . . ,N — 1. A regularized DTH trajectory 
must satisfy the following equations and inequalities: 

Az k = \ k JH z (z k ) + p, k Jip z {z k ) (3a) 

H(z k ) = (3b) 

i>(zk) >i>k (or ip(z k ) < ipk) (3c) 

(J-k (tp(zk) - ipk) = (3d) 

Hk < (or fj, k > 0) (3e) 

Proof. Define the Lagrangian function 

C(zk, Zk+i, Afe, Hk) = A(zk, z k +i) + A fc H(z fc ) + fj,kip(zk) 

where A^ is a Lagrange multiplier for the equality constraint H(zk) = and 
[ik is a Karush-Kuhn- Tucker (KKT) multiplier for the inequality constraint 
ip(zk) > V'fe ( or ip{zk) < ipk)- Applying the KKT necessary conditions, [2], 
jlj, to the regularized DTH principle results in the following equations: 

^Aq k + ^\ k H p (z k ) + -) 
1 . 1. , , 1 



^Pfc = ~o Aa k + ^kH P {zk) + ^fJ-k^p(zk) = (4a) 



Afe+i = 2 Apk + 2 X kHq{Zk) + ^kipq(zk) = (4b) 

£ Afe = H(z k ) = 
V'^fc) > ^k (or V(^fc) < V>fc) 

Hk (V'(^) - ^) = o 

/Ufc < 0, (or \i k > 0). 
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Equations (|4ajl and (|4b|) can be rearranged and combined to give equation 
(EiJ). ■ 

Time reversibility of the regularized DTH trajectory follows from the 
following observation. If the inequality constraint vp(z k ) > ipk is active 
forward in time, the inequality constraint ip(z k ) < V'fc is active backward in 
time. Therefore, the same equation, ip(z k ) = ip k , applies for both directions 
in time. 

3.2 Symplectic-Energy-Momentum Properties 

Next, we show that the regularized DTH equations (|3a|) — (|3e|) is a SEM inte- 
grator. Conservation of energy (Hamiltonian) follows directly from equation 
()3b|) . To prove symplecticness, we identify a generating function which gen- 
erates symplectic transformations between adjacent vertices of a regularized 
DTH trajectory. 

Theorem 7 (Generating Function) Assume [x k and ip(z k ) are not si- 
multaneous equal to zero. Then the function 

S(qk,Pk+i) = qk T Pk+i + £{zk, z k+l, Afc, Hk) 

is a generating function which determines a symplectic transformation be- 
tween the vertices z k and z k+ \ of a regularized DTH trajectory for k = 
0,...,N — 1. The function C(z k , z k +i, ^Jfc> Mfc) is the Lagrangian function 
used in the proof of Theorem . The variables qk+i, Pk, ^k o,nd /j, k are 
determined by the regularized DTH equations \3a\) ~H3e }) . 

Proof. We show that S qk (q k ,p k+1 ) = p k . The equation S Ph+1 (q k ,p k+1 ) = q k 
follows in a similar fashion. 

Sq k = Pk+1 + £q k 

= p k+ i-^Ap k + -XkH q (z k ) + -Hki> q (zk) 

1 A 1 A 

= Pk+l-^pk-^pk 
= Pk 

m 

We note that the transformation z k+ \{z k ) may not be differentiable when 
both if)(z~ k ) and fi k are equal to zero. Theorem may not be valid for this 
special case. 
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The following lemma will be used to prove conservation of momentum. 
(Conservation of momentum is restricted here to linear and quadratic func- 
tions e.g. linear and angular momentum in Cartesian coordinates.) 

Lemma 8 If L{z) is a quadratic function and the Poisson bracket [L,7i] is 
identically equal to zero, then the Poisson bracket [L, ip] is identically equal 
to zero. 

Proof. By assumption, [L,H] = Lj JH Z = J ili2 L h Hi 2 = 0. (We use the 
convention of summing over repeated subscript and superscript indices.) 
Taking the first and second derivative of J lll2 Li 1 Hi 2 = with respect to kth. 
component of z and using the fact that, since L{z) is quadratic, Li lklk2 = 0, 
we have 

J* 1 * 2 Lj 1 fc 1 7Yj 2 + J n * 2 Lj 1 7Yj 2 fc 1 = (5a) 
J 1112 Li 1 k 1 Tii 2 k 2 + J %ri2 Li l k 2 7~Li 2 ki + J %112 Pi- tL Hi 2 k 1 k 2 = 0. (5b) 

Since tp = (JH Z ) T H ZZ {JH Z ) = J kljl 3 k2 i 2 H klk2 H h H h we have 

= J J 2 ^ 2 T~Lk l k 2 i 2 1~Lj l 1~Lj 2 + J 1 ^ 1 J 2 ^ 2 7~L k i k2 7~(-j±i 2 'Hj 2 

+ j k ^j k2 ^ 2 n klk2 n n n j2i2 (6) 

— J J 2 ^ 2 7~ikik 2 i 2 7~tji7~tj2 ^"J J 2 ^ 2 '^-k\k 2 ^-j\i2^j2 

where we combined the last two terms of (jBJ by renaming indices and using 
the fact that Ti k2kl = tik\k%- Substituting for ipi 2 in the Possion bracket 
[L,tp] = Ljjip z = J %1%2 Li^i 2 we have 

M] = j k ^j k2 ^ 2 {j^ 2 L n n klk2i2 ) n n n J2 
+ 2j k ^j k2 > 2 (J hl2 L h n ni2 ) n klk2 n j2 . 

Using equations lj5ajl and (|5b|) and Hj^ = 7~ii 2 j 1 , 7~L klk2 i 2 = r Hi 2 k\ki-i we 
have 

[L, ip] = -J k ^jk2h (j^^H^ + jiM L . ikaHiah ) H n H j2 

- 2J k ^J k ™ (j^ 2 L ilh H l2 ) H klk2 H 32 (7) 



9 



The last term of (JJJ) can be expressed as (jSJ) by rearranging terms, using 
jkiji — _ jjiki anc j renaming indices as shown below. 

_ 2 jk ljljk2j2 (jhi2 L . ijiH .J H klk2 H j2 = 
_ 2J i^jk 2j2 (j"di L . ihHkik ^ Hi2Hj2 = 

2J iii2jk2h (ji^L^H^) H i2 H n = 
2J ^hjk2h (j^ Liiki H l2k2 ) H n H j2 (8) 

Replacing the last term in (JJJ) with (|SJ) and rearranging we have 

[L, iP] = J^hjk^njhn (L hkl H i2k2 - L nk2 H t2kl ) H h H ja (9) 

Skew-symmetry with respect to the indicies k\ and ki in Q implies [L,tp] = 
0. ■ 

Theorem 9 (Quadratic Conservation Laws) Assume L(z) is a quadratic 
function and [L,Ti] is identically equal to zero. Then L(z) is exactly con- 
served at the vertices of a regularized DTH trajectory. 

Proof. Since L(z) is quadratic, L{z k+ i) — L(z k ) = L z (z k ) T Az k . From 
equation Q3a[) and Lemma |H] we have 

L(z k+1 ) - L(z k ) = L z (z k ) T [\ k JH z {z k ) + fi k Jip z (z k )\ 

= {\ k [LM+Hk[LM)\ z= - Zk 
= 0. 



Corollary 10 (Quadratic Conservation Laws) IfL(z) is quadratic, L t and 
L p are both zero, and the Poisson bracket [L, H] is identically equal to zero, 
then L(z) is exactly conserved at the vertices of a regularized DTH trajectory. 

Proof. Since [L,H] = [L,H]+ L t H p — L p TCt, then L t = L p = implies 
[L,7i] = [L,H] and the corollary follows from Theorem |HI ■ 

3.3 Coordinate Invariance 

We briefly consider the coordinate invariance of regularized DTH dynamics. 
In the lemma below, we show that tp(z) is coordinate invariant with respect 
to linear, symplectic, coordinate transformations. 
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Lemma 11 Let z = TZ be a linear, symplectic, coordinate transformation 
between old coordinates z and new coordinates Z. Let Ti(z) be a Hamil- 
tonian function expressed in the old coordinates and IC(Z) = TC(TZ) be 
the Hamiltonian function expressed in the new coordinates. Define ijp^{z) = 
(JH Z ) T H ZZ (JH Z ) and^iZ) = (JJC Z ) T JC ZZ (JJC Z ). Then^{Z) = ^(TZ). 

Proof. Since K{Z) = H(TZ) we have K z = T T H Z and K zz = T T H ZZ T. 

i; K (Z) = (JIC Z ) T IC ZZ (JK Z ) 

= (jt t h z ) t {t t h zz t) (jt t h z ) 
= nj (tj t t t ) n zz (tjt t ) n z 

Since T is symplectic, TJT T = T T JT = J and we have 

iP K (Z) = HJJ T H ZZ JH Z 

= (JH Z ) T H ZZ {JH Z ) 
= i> H {z) 
= ip n (TZ). 



Theorem 12 (Linear Coordinate Invariance) The regularized DTH equa- 
tions are coordinate invariant under linear, symplectic, coordinate transfor- 
mations. 

The proof of Theorem ^] parallels the proof given in [Tl]| for the DTH 
equations. The only difference is the use of Lemma ^2 in the proof of 
Theorem 1121 (We point out that the Lagrange multiplier and the KKT 
multiplier are both coodinate invariant quantities.) 

Ge [0] demonstrated the coordinate invariance of a variety of symplec- 
tic integrators under linear, symplectic coordinate transformations. Gui- 
bout and Bloch demonstrate coordinate invariance using the larger 
class of linear, symplectic, discrete-time coordinate transformations. In fact, 
it should be possible to demonstrate coordinate invariance using the even 
larger class of piecewise-linear, continuous, symplectic coordinate transfor- 
mations which are consistent with a special triangulation of extended phase 
space. A procedure for demonstrating the coordinate invariance of DTH dy- 
namics using this larger class of coordinate transformations was described 
(in a formal sense) in ^2]- Theorem 1121 implies this procedure is also valid 
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for regularized DTH dynamics. In fact, the regularization described in this 
article removes the technical difficulty of identifying principle DTH trajec- 
tories described in [T7] . 

4 Numerical Results 

4.1 An Algorithm for Regularized SEM Integration 

An algorithm for solving the regularized DTH equations (|3a }) -(|5e ]l is outlined 
in Figure In this section, we will choose ipp. = in equation (|3c|) . Before 
explaining the algorithm in detail, it would be useful to review the simpler 
algorithm developed in |17j for solving the DTH equations (|laj) - (|lbj) . 

Equations (|la |) -(|lb [) are poorly conditioned for small time steps A&. A 
direct application of Newton's method is likely to result in poor convergence. 
Instead, nested, Newton iterations are used. The function z^ = z(Xk,Zk) 
implicity defined by equation (jlaj) . is evaluated using an inner iteration. An 
outer iteration solves the equation g(X) = TC(z(X, Zj~)) = for A&. Quadratic 
convergence of the iterations is proved in |17j for V^fc) 7^ 0- The outer 
Newton iteration exhibits poor convergence near ■0 = 0. 

The algorithm outlined in Figure 0] also uses nested, Newton iterations. 
Near tp = however, a bracketed, root-finding procedure is used in place 
of the outer Newton iteration. The algorithm is further complicated by the 
problem of identifying when = has been crossed. The procedural logic 
needed to compute a DTH trajectory crossing ip = appears to be complex. 
We now give a more detailed explanation of the algorithm outlined in Figure 

m 

The function z(A&, fj,f.> z k) implicitly defined by equation (|3"a)l is evaluated 
using Newton's method. When ifi(zk) / 0, equation (|3d|) implies /ifc = 0. 
We use the abbreviation z(\k, z^) for ~z{\k, /Ufc, z^) when ^ = 0. 

The first block in Figure 0] initializes the algorithm. (We assume ip(zo) ^ 
0.) The value of po (the momentum conjugate to time) determines the initial 
time step Ao- If po is chosen so that T~L{zq) = 0, then Ao = 0. Therefore, a 
value for po should be chosen so that TL{zq) is sufficiently small but nonzero. 

Vertex points z^, k = 2, . . . , N, are computed within block 2. To avoid 
ill-conditioning of the equation Tt(z(\, Zf.)) = near ip = 0, we solve 
the equation ^(^(A, Zfc))?Y(z(A, z^)) = 0. Linear segments which do not 
cross ip = are computed in block 3. If ifi(z(\k, Zk)) = 0, the condition 
ip{zk)^p{zk+i) < indicates ■0 = has been crossed and the algorithm en- 
ters block 4 where A^ and z$ are computed and used in block 5. In block 5, 
a bracketed root-finding procedure is used to solve the now ill-conditioned 
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input zo, set k = 

solve ip(z(X, z k ))Tt(z(X, z k )) = for A > to determine X k 
z k = z(X k , z k ), z k+1 = 2z k - z k 
repeat while k < N 

while ip(z k )i/j(z k +i) > and k < N 
k = k + l 

solve i/)(z(\, z k ))H(z(X, z k )) = for A > to determine X k 
z k = z(X k , z k ), z k+ i = 2z k - z k 

end 

solve ip(z(\, z k )) = to determine Ay, 

Zip = z(X^,z k ) 

while Ay, > and H(z k )H(z^) < and k < N 

solve H(z(X, z k )) = for < A < Ay to determine X k 
Zk = z{X k , z k ), z k+1 = 2z k - z k 
k = k + l 

solve ip(z(X, z k )) = to determine Ay, 

z$ = z(X^,z k ) 

end 

if ghost trajectory and H(z k -i)Tt(z k ) > 
k = k- 1 

solve ip(z(X, z k )) = to determine Ay 
solve H(z(X, z k )) = for A > Ay, to determine X k 
z k = z(X k , z k ), z k+ i = 2z k - z k 
k = k+l 

end 

if regularized trajectory 

H(z(X,n,z k )) = 
tp(z(X, n, z k )) = 
z k = z(X k , fi k , z k ), z k+ i = 2z k - z k 
k = k + 1 

solve ip(z(X, z k )) = to determine Ay 
solve TC(z(X, z k )) = for A > max(0, Ay) to determine X k 
z k = z{X k , z k ), z k+1 = 2z k - z k 
k = k + l 

end 

solve H(z(X, z k )) = for A > to determine X k 
Zk = z{X k , z k ), z k+1 = 2z k - z k 
k = k + l 
end 



solve 



to determine X k and fi k 



Figure 4: Algorithm for computing ghost and regularized DTH trajectories. 

13 



x 10 _B x 10" 5 x 10" 




-0.02 -0.01 0.01 0.02 -1 -0.5 0.5 -0.2 0.2 0.4 



(a) Parameter (b) A < 0, /i-small (c) A < 



x 10~ 5 x 10" 5 x 10" 




(d) A = (e) A > (f) A > 0, p-small 

Figure 5: A one parameter family of DTH trajectories crossing ij) = 0. 
Dashed curves are ghost trajectories. Solid curves are regularized trajecto- 
ries. 

equation 7i(z(X, Zk)) = 0. If a bracket can not be found, the algorithm 
enters either block 6 and computes a ghost trajectory or block 7 and com- 
putes a regularized trajectory. Block 8 is need to prevent the computation 
of trajectories which immediately recross ip = 0. 

Blocks 7 and 8 take into account the different ways DTH trajectories can 
cross t/j = 0. These blocks are best understood after viewing an animation 
of a one-parameter family of ip = crossings. Snapshots of this animation 
are given in Figure |SJ 

Murua JB] has developed an efficient iteration which avoids the nested 
iterations use to solve 7i(z(X, z^}) = 0. (Murua has also developed an itera- 
tion which does not require evaluation of the Hessian matrix of the Hamil- 
tonian function.) It is likely that the algorithm outlined in Figure 0] could 
be made significantly more efficient by using Murua's iteration. The author 
has take a first step in this direction by modifying Murua's iteration so that 
it can be used to compute regularized segments crossing ip = 0. 
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(a)g = 0,p = 3 (b)g = 0,p = 2.2 

Figure 6: Behavior of A vs /x when -0 = for a single regularized DTH 
trjectory of the nonlinear pendulum. 

4.2 Qualitative Behavior of Regularized SEM Integation 

Numerical computations for the pendulum and Kepler's one body prob- 
lem confirm that the regularization described in this article conserves SEM 
properties. The energy (Hamiltonian) is conserved to roundoff error at mid- 
points of DTH trajectories. Angular momentum is conserved to roundoff 
error at vertices for Kepler's problem in Cartesian coordinates. Symplec- 
ticness is verified by computing the derivative dzjy/dzo of the map zn(zq). 
The matrix (cIzn /dzo) T J (dzjy / dzo) is found to equal J to roundoff error. 
Time-reversibility is also confirmed to hold to roundoff error. Numerical 
computations quantifying the accuracy and efficiency of the regularization 
have not yet been completed. 

One of the peculiarities observed in regularized SEM integration is the 
occurrence of negative time steps. Negative time steps violate the monotonic- 
increasing property of time. The DTH trajectories become multi-valued 
functions of time. Lee ^H] foresaw this possibility and suggested the rem- 
edy of relinking vertices to maintain the monotonicity of time. 

Finally, another peculiarity observed for regularized SEM integration, 
for the case ipk = in equation (jHcj) . is apparently chaotic behavior near the 
separatrix of the pendulum. (See Figure E3) It may be possible to regularize 
DTH dynamics further by choosing nonzero values for -0^. 
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